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Abstract-We propose a method of simulation that is based on the averaging of formal solutions of the transfer 
equation by taking the integral by the Monte Carlo method. This method is used to compute two models, which 
correspond to the limiting cases of hot gas and cold background radiation and of hot background radiation and 
cold gas, for E-methanol emission from a compact homogeneous spherical region. We analyse model level pop- 
ulations by using rotational diagrams in the limiting cases mentioned above. Model optical depths of the lines 
with frequencies below 300 GHz up to J=ll inclusive are given. 



INTRODUCTION 

Molecular methanol (CH3OII) emission, along with 
emission of other molecules, is commonly observed 
toward star-forming regions. Some molecular methanol 
transitions produce narrow and intense maser lines in 
some sources, while, in other sources, these transitions 
produce less intense or no detectable lines; at the same 
time, maser activity may show up in other transitions. 
The conditions for the formation of methanol masers 
can be studied by numerically simulating radiative 
transfer. In such problems, we simultaneously solve 
the transfer equation and the system of statistical- 
equilibrium equations for a large number of lines. This 
problem is commonly solved either by the Monte Carlo 
(MC) method or by large velocity gradient (LVG) 
method. In contrast to the LVG method, the MC 
method allows computations to be performed for a 
compact cloud and a small velocity gradient. In ad- 
dition, the model adopted in the MC method admits 
several important complications: inhomogeneity and 
complex structure of an emitting cloud and a more 
realistic allowance for the effect of its nonspherical 
geometry (with no use of factor e^^). However, the MC 
method requires more computing time than the LVG 
method does. The MC method, whose algorithm was 
described by Bernes (1979), is based on the replace- 
ment of the real radiation field by a number of model 
photons and on the simulation of their propagation 
through the medium with the computation of the num- 
ber of molecular excitations in each of the shells into 
which the cloud is divided. The classical MC method 
assumes that model photons are emitted uniformly in 



all directions in each of the shells and come from 
the outside. Juvela (1997) proposed a modification of 
the standard MC method: the paths of all model pho- 
tons begin at the cloud edge. However, as in the classical 
MC method, the formalism of model photons is used as 
the basis for the computations. In this paper, we make 
an attempt to deduce the method from the transfer 
equation while keeping the main features of the MC 
method and propose a method that is based on the 
averaging of formal solutions to the transfer equation by 
taking the integral by the Monte Carlo method. This 
method was tested on simulations of E-methanol emis- 
sion from a compact (~0.005 pc) spherical cloud with 
no velocity gradient. 

THE ALGORITHM AND BASIC 
APPROXIMATIONS 



We performed simulations in the simplest case of a 
static spherical cloud with a uniform distribution of the 
H2 and CH3OH number densities and kinetic tempera- 
ture. The influence of dust was ignored. The spectrum 
of the background radiation was assumed to follow the 
Plank law. The dilution of the background radiation 
was disregarded. The radius of the cloud was set equal 
to 1.5x10^6 cm. 

In the algorithm under consideration, as in the algo- 
rithm of Bernes (1979), the solution of the problem 
reduces to an iterative procedure which is implemented 
using the system of statistical-equilibrium equations 
(|^) and which yields new approximations for the popu- 
lations 
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where Skj = 1 when and only when the relation between 
the energy levels Ek > Ej holds; otherwise, 6kj = 0; Ajk 
and Bjk are the Einstein coefficients (which are assumed 
to be zero for forbidden transitions); Cjk are the colli- 
sional constants; nk is the population of the kth level; 
and / is the average intensity. In the first equation in 
(|I]), the subscript j changes in the range 1 < j < (A^ — 1), 
where N is the total number of levels involved. The 
average intensity is given by 



+ 00 



+1 



i{v, ^) d^i, 



(2) 



where v is the frequency, is the cosine of the angle 
between the radius and the selected direction, and f[v) 
is the line profile (here it was assumed to be the Dop- 
pler profile). This integral can be calculated by the 
Monte Carlo method by simulating the random vari- 
ables u and /i with Gaussian (for the Doppler line pro- 
file) and uniform distributions, respectively: 



— -l = - fcrsin(27ri?) {-\nR'f'^+v[r)ii 
Vq c \ 



^i = ^- 2R" 



(3) 
(4) 



where R, R' and R" are independent random variables 
with a uniform distribution in the interval to 1; z^o is 
the rest frequency of the transition under considera- 
tion; v{r) is the velocity field at a given point; and a 
is the Doppler halfwidth. The average intensity can be 
replaced by the mean 



/ = MI {fi, ly) 



(5) 



The averaging over sufficiently thin shells into which 
the cloud is divided leads to an additional integration of 
(H) over the radius, which reduces to the simulation 
of yet another random variable that gives the distance 
from the cloud center to the common point of the bun- 
dle of directions used to compute the mean (H) : 



1/3 



(6) 



where rj„ and rout are the radii of the inner and outer 
boundaries of the shell, and i? is a random variable that 
is uniformly distributed in the interval to 1 . The inten- 
sity I {fijv) for a given direction can be estimated by 
using 



L (m) = L bg {^J■) exp I - / x^{y) dyj 
+ f £^(X) exp i - / Xu{y) dy \ dX, 





(7) 



where and Ei, are the absorption and emission coeffi- 
cients of the medium, respectively; L is the distance to 
the cloud edge in the direction of integration; and bg is 
the intensity of the background emission for a given 
direction. The integration path is broken up into seg- 
ments in which the emission and absorption coeffi- 
cients can be assumed to be constant (in the presence of 
a velocity gradient, this path can be much shorter than 
the path within a single shell where the populations are 
assumed to be constant). The change in intensity within 
a segment can be obtained by integrating (^: 



1 



(8) 



where is the intensity at the beginning of the seg- 
ment, and I is the segment length in the direction under 
consideration. The emission and absorption coeffi- 
cients are expressed in terms of the level populations 
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(10) 



where is the Einstein coefficient for the correspond- 
ing transition; n„, and n/, gi are the populations and 
statistical weights of the upper and lower levels, re- 
spectively. Thus, the average intensity is expressed in 
terms of the level populations and the intensity of the 
background radiation, which allows us to implement the 
iterative process using the system of statistical-equilib- 
rium equations (0) for the populations. When the itera- 
tive process has converged, the derived populations are 
used to compute the intensity of the cloud radiation (or 
brightness temperature) with the aid of the formal solu- 
tion to the transfer equation (0). 

RESULTS 

We performed the computations for two models 
with kinetic temperatures of 70 and 20 K and with 
background-radiation temperatures of 2.7 and 70 K, 
respectively (below referred to as models I and II 
respectively). The methanol column densities in the 
two models were assumed to be the same and equal to 
1.5 X 10^^ cm~^ (the radius is 1.5 x 10^^ cm, the II2 num- 
ber density is 10^ cm"-^, and the methanol abundance is 
10~^). When computing the models, we took into 
account 124 lower rotational levels of the ground tor- 
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sional state of E-methanol (J < 15, \K\ < 6, and 
E < 200 cm-i). 

Torsional transitions were ignored. We computed 
the energies of the levels and Einstein coefficients 
(A) by using approximate formulas from Pikett et al. 
(1981). The collisional constants are known poorly and 
were computed here in the gas-dynamical approxima- 
tion of Lees (1974). The iterations were terminated at 
an accuracy of 0.1% of the population. 

DISCUSSION 

The algorithm of Bernes (1979) uses the quantity 
Siu,m' (for the m' shell and for the transition between / 
and u) which is proportional to the intensity instead of 
the latter; this quantity is accumulated over all model 
photons in all shells. A comparison of the system of 
statistical-equilibrium equations yields the relation 
(X) Siu,m')—Iiu,m'Biu- In each step, according to formula 
(6) from Bernes (1979), the following quantity is added 

to (S Slu,m')- 



Slu,r. 



SkWo exp I - 

= ^f{^)Biu yj.-^ ' {1 - cxp (-r,.)}. 



(11) 



where Wq is the initial weight of the model photon; Vm' is 
the volume of the m' shell; and and Si are the optical 
depth and path of the model photon in the «th step, 
respectively. In order to compare the method of Bernes 
(1979) with the method proposed here, which we 
deduced from the transfer equation, let us consider the 
contribution of the intrinsic emission from the m shell, 
the emission from a different m' shell, and the back- 
ground radiation to (^ Siu,m')- Substituting the relations 
for Ti and for the initial weight of the model photon 
produced in the m' shell, WQ=Vm' A^iUu / Nph, where Nph 
is the total number of model photons, we obtain the fol- 
lowing expression for the contribution of the intrinsic 
emission from the m' shell: 



c, _ Biu 



{l-exp(-Ti)}. (12) 



Since the term in parentheses is the source function and 
comparing this equation with (H), we conclude that the 
contribution of the intrinsic emission is accurately rep- 
resented in the method of Bernes (1979). It follows 
from (|^) that the contribution of the background radia- 
tion must be given by 



fe-i 



Slu,ra' = -j^^bg CXp - - 



(13) 



where Ihg is the intensity of the background radiation at 
the frequency of the transition between levels u and 
I. If the cloud is divided into shells of equal volume 
(Kn' — const) ^ then the path of the model photon is 
the same in each step (sfc — const) and small enough 
for the optical depth in each step to be small in ab- 
solute value. Formula (|l^) can then be derived from 
( pj] ) by using the linear approximation for the exponent. 
The contribution of the emission from a different shell 
can be obtained from ( pi] ) , by analogy with ( p^ ) , 



W„h v„, '-^V 



fc-i 

En 

i=l 



—pr- 



{1 -CXp(-Tfe)}, 



(14) 



where the subscript k of the source function means that 
this function is calculated from the populations in the 
fcth step (in the m' shell rather than the m shell where 
the model photon was formed). It follows from (^ that 
this contribution must be equal to the product of (|l2|), 
where the source function is calculated from the popu- 
lations in the first step and the exponential factor from 
(p^). If the cloud is divided into shells of equal volume 
and if the optical depth in each step is small in abso- 
lute value (which allows a linear expansion of the expo- 
nent), then, using relation (^) for the emission coeffi- 
cient, we obtain the last condition for the applicability 
of the algorithm of Bernes (1979): 



fk{v){nu)kSk ~ fi{y){nu)isi, 



(15) 



where the subscripts denote the first and fcth steps in the 
propagation of the model photon. Relation ( p^ is the 
condition for the equality of the specific column densi- 
ties of the emitting molecule in each step in the direc- 
tion of propagation of the model photon. The latter con- 
dition is most stringent and is difficult to satisfy in prac- 
tice. Since uncertainty in the estimate of the radiation 
field in the algorithm of Bernes (1979) produces addi- 
tional noise in the method, the use of the algorithm out- 
lined here seems more appropriate. Our method cor- 
rectly describes the radiation field (the formulas follow 
from the transfer equation); therefore, its application 
is restricted only by the convergence of the iterative 
procedure used (matches the iterative procedure in the 
classical MC method) and by the available computing 
time. In practice, the convergence can be hampered in 
the presence of strong masers. 

The qualitative behavior of populations in the limit- 
ing cases of hot gas (model I) and hot background radi- 
ation (model II) is convenient to analyze as follows. Let 
us consider the relation between the ratio of the level 
population to its statistical weight and the energy of this 
level on a logarithmic scale. In order not to overload the 
picture, let us consider only the ladders with quan- 
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Fig. 1. Logarithm of the ratio of the population (in cm~^) to 
the statistical weight versus the energy of a given level for 
the model with a background-radiation temperature of 2.7 K 
and a cloud kinetic temperature of 70 K. The dashed lines 
correspond to rotational temperatures of 64 and 28 K. 

turn numbers K = and -1. In the LTE case, the 
level population would be given by Boltzmann formula, 
and we would have a straight line whose slope would be 
related to the temperature. In the non-LTE case, how- 
ever, we have a curve (Figs. 1 and 2). The levels that 
correspond to the different ladders in these figures are 
indicated by different symbols. In the model with hot 
gas, the K = —1 ladder is heavily overpopulated relative 
to the K = ladder at small quantum numbers J, while 
in the radiation-dominated model, the K = ladder lies 
above the K = —1 ladder; i.e., there is an inversion of 
other transitions. This may correspond to the division 
of methanol masers into class I (Fig. 1) and class II 
(Fig. 2) masers (Batrla et al. 1987; Menten et al. 1991). 
The slope of the curve at energies E < 100 cm~^ in Fig. 1 
corresponds to a rotational temperature T^ot 64 K, 
which is close to the kinetic temperature of the 
medium. At high energies, T^ot ~ 28 K. In Fig. 2, rota- 
tional temperatures Trot ~ 27 K (which is close to the 
kinetic temperature of the medium) and Trot ~ 59 K 
(which is close to the background-radiation tempera- 
ture Tbg = 70 K) correspond to low and high energies, 
respectively. Thus, the slope of the curve appears to be 
mainly determined by collisions at low energies and by 
radiation at high energies. The distribution of popula- 
tions within the ladder (in each of the segments) is close 
to the Boltzmann distribution, while, between the lad- 
ders, it is an purely non-LTE distribution. For this rea- 
son, no emission arises in a- type transitions {AK = 0). 



Fig. 2. Same as Fig. 1 for the model with a background-radi- 
ation temperature of 70 K and a cloud kinetic temperature 
of 20 K. The dashed lines correspond to rotational tempera- 
tures of 27 and 59 K. 

The b-type transitions (which occur with a change in 
the quantum number K) tend to produce a series of the 
form {J+a)K±i—JK, where a = 0, ±1 (Menten et al. 
1986b). Since the scries Jq — (J+l)-i includes transi- 
tions from the K=0 ladder to the K=—l ladder and vice 
versa, some of the transitions in this series exhibit activ- 
ity in model I (4_i — 3o, 5_i — 4o, etc.), and other arc 
active in model II (Oq — 1-1, Iq — 2_i, 2o — 3_i). Masers in 
this series were detected in the observations of Turner 
et al. (1972), Zuckerman et al. (1972), Batrla et al. 
(1987), and Slysh et al. (1997, 1999). In the models 
under consideration, the series Jq— J-i at 157 GHz gives 
a weak inversion of transitions for J=2 and 3. Masers 
in this series were detected in the observations for 2< J<8 
(Slysh et al. 1995). The transitions in the series Jq— 
(J-l)_i have higher frequencies and show no inversion 
in the models under consideration. A similar behavior 
of the populations is observed for transitions between 
other ladders. The (J+l)o— Ji transitions exhibit class I 
and II activity for J > 2 and J < 3, respectively. Wilson 
et al. (1985) detected a maser in the 2i — 3o line at 
20 GHz. Slysh et al. (1992) identified weak thermal 
emission in the 4o — 3i line at 28 GHz. The series Ji— Jo 
at 165 GHz and Ji — (J-l)o give no inversion in the 
models under discussion. No masers in the series Ji— Jq 
were detected in the observations of Slysh et al. (1999). 
The transitions in the series J2— Ji at 25 GHz give an 
inversion in model I. No masers in this series were 
detected in the observations of Barret et al. (1971) and 
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Peak optical depths of the hnes for models with the following H2 number densities, methanol abundances, kinetic 
temperatures, and background-radiation temperatures: 1 - = 10^ cm~^, X = 10^, T^m = 70 K, and Tf,g = 2.7 K; 
2-nH2= 10^ cm-3, X = 10^ Tki„ = 20 K, and Tbg = 70 K 
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25.9 


-0.7* 




0.0 


7o - 7-1 


156.8 


7.7 


M2'' 


0.7 


11-5 - 12-4 


172.1 


0.0 




0.0* 


7i - 7o 


166.2 


3.7 




0.8 


11-3 - 10-4 


7.3 


0.0 




0.1 


72 - 7i 


25.1 


— 1.6* 


Ml*^ 


0.7 


11-3 - 12-2 


183.7 


0.1 




0.0 


'3—02 




i.i 




— (J.i 


11 in 

11-1 — 10-2 


1 H/l Q 

iU4.o 


— u.y 




v. 6 


8-4-9-3 


89.5 


0.0 




-0.1* 


llo-ll-i 


154.4 


2.4 




0.3 


8-1 - 7o 


229.8 


-1.0* 


Ml" 


2.0 


111 - llo 


171.2 


0.8 




0.3 


80 -7i 


220.1 


-0.1* 




1.4 


112 - 103 


2.9 


-0.3* 




0.1 


80 — 8-1 


156.5 


6.5 


M2'' 


0.6 


112-111 


26.3 


-0.4* 




0.0* 


81 — 80 


166.9 


3.2 




0.6 


114 - 123 


253.9 


0.0 




0.0 



Note; Ml - a class I mascr was observed; M2 — a class II mascr was observed; — M - no mascr was detected; * - the optical depth in this 
transition is negative; (a) - Slysh et al. (1999); (b) - Slysh et al. (1995); (c) - Batrla ei al. (1987); (d) - Wilson et al. (1985); (e) - 
Barret et al. (1971); (f) - Turner et al. (1972); (g) - Slysh et al. (1992); (h) - Slysh et al. (1997); (i) - Zuckerman et al. (1972); (j) - 
Haschick et al. (1989); (k) - Slysh et al. (1993); (1) - Mcntcn et al. (1986a). 



Menten et al. (1986a). Inversion in class I models is also 
observed for the transitions of the series (J+l)2— Ji- 
The first two transitions in the scries ( J = 1 and 2) turn 
out to show the greatest inversion. In model 11, inver- 
sion is observed in the 3i — 22 transition from the series 
(J+l)i— J2. Since the series J-2 — (J+l)-i contains 
transitions from K = —2 ladder to the K = —1 ladder 
and vice versa, some of the transitions show inversion 
in model I (9-i — 8-2. lO-i —9-2, etc.), and others show 
inversion in model II (5-2—6-1, 6-2—7-1, 7-2— 8-1; tran- 
sitions with smaller J have higher frequencies and show 
no inversion or their inversion is marginal). The lowest 
frequency transitions 7-2 — 8-1 at 37 GHz and 9-1—8-2 
at 9.9 GHz, which were observed in the direction of 
star-forming regions (Haschick et al. 1989; Slysh et al. 
1993), are the most intense transitions in this series. No 



maser emission was detected in the 5-2—6-1 and 3-2—4-1 
transitions (Slysh et al. 1999). The series (J-|-l)2— J3 
(J> 9) gives inversion in model I. The optical depth in 
transitions between ladders with different K is small 
because of the low population (although inversion may 
exist). Our computations arc in qualitative agreement 
with the data of Cragg et al. (1992), which were 
obtained by the LVG method. Candidates for class II 
mascrs coincide with those of Sobolev et al. (1997), 
with the exception of the series Ji— Jo at 165 GHz. This 
result may be related to the difference in the model 
parameters. Peak optical depth for the lines with fre- 
quencies below 300 GHz up to J = 11 inclusive in models 
I and II are given in the table. When computing the 
spectra, we assumed that all sources were observed 
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against the background radiation with a temperature of 
2.7 K. 

CONCLUSION 

(1) The method of Bernes (1979) is apphcable if the 
cloud is divided into sheUs of equal volume, if the paths 
of model photons and the optical depths in each step are 
sufficiently small, and if the column densities of the 
emitting molecule in the direction of propagation of the 
model photon are equal for each step. The applicability 
of the method proposed here is restricted only by the 
convergence of the iterative procedure, which may be 
hampered in the presence of strong masers. 

(2) The populations within the ladder (Figs. 1 and 2) 
in the segments before and after the break are well 
described by the Boltzmann formula, while the popula- 
tions between the ladders have a non-LTE distribution. 
For this reason, no masers are formed in transitions that 
occur with no change in the quantum number K. The 
position of the break in the "logarithm of population- 
to-statistical-weight ratio versus level energy" diagram 
appears to be determined by the ratio of the rates of col- 
lisions and radiative processes. 

(3) The transitions 22 - li at 121 GHz, 82 - 2i at 
170 GHz, IO2 - lOi and II2 - Hi at 26 GHz, 5o - 4i at 
76 GHz, 60 - 5i at 124 GHz, 7o - 61 at 172 GHz, 7_i - 60 
at 181 GHz, 10_i - 9_2 at 57 GHz, ll_i - 10_2 at 
104.3 GHz, II2 - IO3 at 2.9 GHz are candidates for 
new class I masers, while the transitions Iq — 2_i at 
61 GHz, li - 2o at 68 GHz, and 6_2 - l-i at 86 GHz are 
candidates for class II masers. 
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